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Abstract 



By employing a semi-analytical dynamical mean-field approximation theory 
previously proposed by the author [H. Hasegawa, Phys. Rev. E 67, 041903 
(2003)], we have developed an augmented moment method (AMM) in order 
to discuss dynamics of an iV-unit ensemble described by Langevin equations 
with delays. In AMM, original A-dimensional stochastic delay differential 
equations (SDDEs) are transformed to infinite-dimensional deterministic DEs 
for means and correlations of local as well as global variables. Infinite-order 
DEs arising from the non-Markovian property of SDDE, are terminated at the 
finite level m in the level-m AMM (AMMm), which yields (3+m)-dimensional 
deterministic DEs. Model calculations have been made for linear and nonlin- 
ear Langevin models. The stationary solution of AMM for the linear Langevin 
model with N = 1 is nicely compared to the exact result. In the nonlinear 
Langevin ensemble, the synchronization is shown to be enhanced near the 
transition point between the oscillating and non-oscillating states. Results 
calculated by AMMO are in good agreement with those obtained by direct 
simulations. 
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I. INTRODUCTION 



The time delay plays an important role in many systems such as optical devices [1], 
physiological [2] and biological systems [3] . The effect of time delays has been theoretically 
studied by using the time-delay differential equations (DEs). Its exposed behavior includes 
the multistability and the bifurcation leading to chaos. It is well known that noise also plays 
important roles in these systems, and its effects have been thoroughly investigated with the 
use of stochastic DEs. One of its representative phenomena is the stochastic resonance [4], 
in which the signal-to-noise ratio is enhanced for sub-threshold signals. 

In real systems, both noises and time delays coexist, and the combined effect may be 
described by stochastic delay differential equation (SDDE). For instance, SDDEs are used 
in optics [5] and physiology [6] to model noise-driven systems exhibiting delay feedback. 
In recent years, there has been a growing interest in combined effects of noise and delay. 
The theory for SDDE remains much less studied and has been a subject of several recent 
papers [7]- [14], in which the stability condition for the equilibrium solution of linear delay 
Langevin equation has been studied. Its stationary solution is investigated by using the step 
by step method [7] and the moment method [8] . The Fokker-Planck equation (FPE) method 
is applied to SDDE in the limit of a small delay [9]. These studies have been confined to 
the stationary solution of SDDE. More interesting is, however, expected to be its dynamics 
in a stochastic system with a large time delay. 

Real physiological and biological systems usually consist of many elements, each of which 
is described by SDDE. A typical example is a living brain, in which a small cluster contains 
thousands of similar neurons. Each neuron which is subject to various kinds of noises, 
receives spikes from hundreds of other neurons through dendrites with a transmission delay 
and generates spikes propagating along axons. Theoretical study on such coupled, stochastic 
systems has been made by using direct simulations (DSs) [15-17] and analytical methods 
like FPE [18]. Since the time to simulate such systems by conventional methods grows as iV 2 
with N, the size of the ensemble, it is rather difficult to simulate systems with the realistic 
size of N ~ 100 — 1000. Although FPE is a powerful method in dealing with the stochastic 
DE, a simple application of FPE to SDDE fails because of the non-Markovian property of 
SDDE: an evaluation of the probability density at the time t requires prior knowledge of the 
conditional probability density between times of t and t — r, r being the delay time. 

Quite recently the present author [19] [20] has proposed a dynamical mean-field approx- 
imation (DMA) as a semi-analytical method dealing with large-scale ensembles subject to 
noises, extending the moment method [21]- [23]. DMA has been first applied to an N- 
unit ensemble described by the FitzHugh-Nagumo (FN) neuron model without time delays 
[24] , for which original 2iV-dimensional stochastic DEs are transformed to eight-dimensional 
deterministic DEs for moments of local and global variables [19]. In a subsequent paper 
[20], DMA is applied to an iV-unit general neuron ensemble, each of which is described by 
coupled ^-dimensional DEs, transforming ifiV-dimensional DEs to A^-dimensional DEs 
where N eq = K(K + 2). In the case of the Ho dgkin- Huxley (HH) model with K — 4 [25], 
we get N eq = 24. The spiking-time precision and the synchronization in FN and HH neuron 
ensembles have been studied as functions of the noise intensity, the coupling strength and 
the ensemble size. The feasibility of DMA has been demonstrated in Refs. [19] and [20]. 

The purpose of the present paper is to apply DMA to Langevin ensembles with delays, 
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which are expected to be good models representing not only interconnected neural networks 
but also social and technological ones. When DMA is applied to ensembles described by 
linear and nonlinear Langevin equations with delays, the original iV-dimensional stochastic 
DEs are transformed to the infinite-dimensional deterministic DEs for means and correlation 
functions of local and global variables. Infinite-order recursive DEs arising from the non- 
Markovian property of SDDE, are terminated at the finite level m in our approximate 
method, which is hereafter referred to as the augmented moment method (AMM). We may 
study dynamics and synchronization of linear and nonlinear Langevin ensembles with delayed 
couplings, and examine the validity of AMM whose results are compared to results of DSs. 
In particular, for the linear Langevin model with N — 1, a comparison is possible with the 
exact stationary solution [7]. 

The paper is organized as follows. In the next Sec. II, we describe the adopted model and 
method to derive the infinite-dimensional deterministic DEs from the original V-dimensional 
stochastic DEs. Infinite-order recursive DEs are terminated at the finite level m in the level- 
m AMM (AMMm). Model numerical calculations for the linear Langevin model are reported 
in Sec. IIIA, where calculated results of AMM6 for N = 1 are nicely compared to exact 
solutions available for the stationary state [7]. Our AMM is compared also to the small-delay 
approximation (SDA) [9] which is valid for a very small delay. In Sec. IIIB, we present model 
calculations for the nonlinear Langevin model in which the stable oscillation is induced by an 
applied spike for an appreciable delay. The synchronization in the ensemble is investigated. 
It is shown that results of AMM6 are in good agreement with those of DSs. The final Sec. 
IV is devoted to conclusions and discussions. In a following paper [26], our AMM has been 
applied to ensembles described by the noisy FN neuron model with delayed couplings. 

II. ENSEMBLES DESCRIBED BY LANGEVIN MODEL 
A. Basic formulation 

Dynamics of a Lanvegin ensemble with delayed couplings is assumed to be described by 

^jp = FMt)) + ^ E H Mt - r)) + m + / (e) (t), (< = 1 to N) (1) 

with 

&\t) = AG(t - t m )Q(t m + T W - t). (2) 

Here F(x) and H(x) are functions of x, whose explicit forms will be shown later [Eqs. (18), 
(19), (29), (30), (43) and (44)]. We have assumed uniform all-to-all couplings of w and time 
delays of r. The former assumption has been widely employed in many theoretical studies. 
The latter assumption may be justified in certain neural networks [27]. White noises of 
are given by < >= and < >= P 2 SijS(t — t') with the noise intensity of f3 

[28]. An applied input of F e >(t) given by Eq. (2) triggers oscillations in ensembles when 
model parameters are appropriate as will be shown in Sec III, Q(t) denoting the Heaviside 
function, A the magnitude, t in the input time, and T w the spike width. 
In DMA [19], the global variable is given by 
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*(*) = ^E**(*), (3) 

i 

with which we define means and correlation functions given by 

p(t) = < X(t) >, (4) 
7(M')4S<^)M0>> (5) 

i 

p(^t') = < SX(t) SX(t') >, (6) 

using 5xi(t) = Xi(t) — p{t) and 5X(t) = X(t) — p{t). 

In deriving equations of motion of means and variances, we have assumed that the 
noise intensity is weak and that the state variables obey the Gaussian distributions around 
their means, as in Refs. [19,20]. Numerical simulations have shown that for weak noises, 
the distribution of the state variable of an active rotator model nearly obeys the Gaussian 
distribution, although for strong noises, its distribution deviates from the Gaussian [23]. 
Similar behavior has been reported also in FN [22] [23] and HH neuron models [29] [30]. 

After some manipulations, we get DEs for p(t), ^(t,t) and p(t,t') given by (for details 
see Appendix A), 

' hi[!) -go(t)+wuo(t-T) + lV(t), (7) 



dt 
dl(t,t) 

dt 



2^(t) 7 (t, t) + 2w Ul (t - r)p(t, t-r)+ [3 2 , (8) 



= 2 9l (t)p(t, t) + 2w Ul (t - r)p(t, t-r) + ^, (9) 

dp(t, t — mr 



dt 



with 



= [gi(t) + gi(t — mr)]p(t, t — mr) + wu\{t — (m + l)r)p(t, t — (m + l)r) 
+ wu±(t — r)p(t — r,t — mr) + — A(mr), (for m > 1) (10) 



n=0 n - \ Z / 



oo 



n=0 U - 



oo 



n=0 nl 



oo 



„ j) = E ^f2MV, (14 , 



n=0 n! 



where A(x) = 1 for x = and otherwise. Equations (7)-(10) show that an equation 
of motion of p(t,t — r) includes new terms of p(t — T,t — r) and p(t, t — 2r), which arise 
from the non-Markovian property of SDDE. The recursive structure of DEs for p(t, t) is 
schematically expressed in Fig. 1, where arrows express the mutual dependence: p(t,t) 
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depends on p(t, t — r), and p(t, t — r) depends on p(t, t — 2r) and p(£ — r, i — r), and so on. 
Then DMA transforms the original A-dimensional SDDEs given by Eqs. (1) and (2) to the 
infinite-dimensional deterministic DEs given by Eqs. (7)- (10). 

In actual numerical calculations, we will adopt the level-m approximation (AMMm) in 
which DEs are terminated at the finite level m: 



p(t,t- (m+l)r) = p(t,t-mr). (15) 

As will be shown later in model calculations with changing m, the calculated result converges 
at a rather small m [Figs. 2(a) and 9(b)]. 

We note that the noise contribution is f3 2 in Eq. (8) and it is j3 2 /N in Eq. (9). It is easy 
to get 

p(M) = ^M (16) 
= 7 (M). hr(3 2 /w-^0 (17) 

Equation (16) is consistent with the central-limit theorem. We will show later that with 
varying model parameters, the ratio of p(t,t)/~f(t,t) is varied, which leads to a change in 
the synchronization of ensembles [Eq. (27)]. 

DSs have been performed for Eqs. (1) and (2) by using the fourth-order Runge-Kutta 
method with a time step of 0.01. Initial values of variables at t G (— r, 0] are Xi(t) = x* 
for % — 1 to N, where x* is the stationary solution for f3 — 0. The trial number of DSs 
to be reported in the next Sec. Ill is N r = 100 otherwise noticed. AMM calculations 
have been performed for Eqs. (7)-(10) with Eq. (15) by using also the fourth-order Runge- 
Kutta method with a time step of 0.01. Initial values are p,(t) = x* and 7(M) = at 
t e [-r, 0], and p(t,t') = at t G [-r, 0] and t' G [-r, 0] (t > t'). All calculated quantities 
are dimensionless. 



III. MODEL CALCULATIONS 



A. Linear model 



We first consider the linear (L) model given by 

F(x) = -ax, (a > 0) (18) 

H(x)=x. (19) 

The stability of the stationary solution of Eq. (1), (18) and (19) with N = 1 and I^(t) = 
was discussed in Refs. [7]- [13], in particular, with the use of the moment method by Mackey 
and Nechaeva [8]. When Eqs. (18) and (19) are adopted, Eqs. (7)-(10) become 

^fi = _ a/i ( t ) + Wfl ( t _ T ) + /(e) ( t ) > (20) 

= -2 ai (t,t)+2wp(t,t-T)+(3 2 , (21) 



dt 
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= _ 2ap (t, t) + 2wp(t, t-r) + £ (22) 

dp(t, t — tot) 



dt 



= —2ap(t, t — mr) + wp(t, t — (to + l)r) 
+ iup(t - r, t - tot) + | y J A(tot), (for to > 1) (23) 



because g (t) = —ap(t), gi(t) = —a, uo(t) = pit) and ui(t) = 1 in Eqs. (11)-(14). 

For P — and I^ e \t) = 0, Eq. (20) has the stationary solution of p* = 0. Linearizing 
pit) around p*, we get the condition for the stationary solution given by 

cos^ia/w) 

r < Tc = » (24) 

which is just the same as the N = 1 case [8]. 
N=l case 

First we discuss model calculations for N — 1, for which the exact solution of its sta- 
tionary state is available. From Eqs. (21) and (22), we get p(t,t f ) = "f(t,t') in the case of 
N — 1. Solid curves in Fig. 2(a) express 7*, the stationary value of j(t,t), when the level 
to in AMMm is varied for a — 1, w — 0.5 and (3 = 0.001, whereas dashed curves denote the 
exact result given by [7] 

„ / it; sinh(rd) - d \ 2 
\ 2d[w coshird) — a\J 



where d = Va 2 - w 2 . Equation (25) yields 10 6 7* = 1-0, 0.724, 0.635, 0.581 and 0.577 for 
r = 0, 1, 2, 5 and 10, respectively. Figure 2(a) shows that for r = 0, the result of AMM 
agrees with the exact one for to > 0. In the case of r = 1, the result of AMM is larger than 
the exact one for to = 0, but the former is smaller than the latter for to > I. This is the 
case also for r = 2. In contrast, in cases of r = 5 and 10, the results of AMM are in good 
agreement with the exact ones for to > 1. It is surprising that the results of AMM converge 
at a small to (~ 1). Solid and dashed curves in Fig. 2(b) show the r dependence of 7* of 
AMM6 and the exact result, respectively (hereafter we show results in AMM6). The result 
of AMM is in a fairly good agreement with the exact one for r > 4. 

Figure 2(b) shows the r dependence of 7*. The result of DMA is in good agreement with 
the exact one for r > 3. It is interesting to make a comparison with results calculated by 
the small-delay approximation (SDA) initiated in Ref. [9], some details of SDA being given 
in appendix B. The dotted curve in Fig. 2(b) expresses 7* calculated in SDA for w = 0.5, 
(3 = 0.001 and N — 1. Although the result of SDA agrees with the exact one at very small 
t (~ 0), it shows a significant deviation from the exact one at r > 2 where 7* becomes 
negative violating its positive definiteness. 

Solid and dashed curves in Fig. 2(c) express the w dependence of 7* of AMM and 
exact ones, respectively, for various r values with (5 = 0.001 and N — 1. We note that an 
agreement between AMM and exact results is good except for w > 0.6 with r = 1 ~ 2 and 
for w < —0.8 with r = 2 ~ 10. From the results shown in Figs. 2(a)-2(c), we may say that 
AMM is a good approximation for a large r (> 4) and a small (3. 
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The response of the Langevin model (N = 1) will be discussed to an applied spike of 
&\t) given by Eq. (2) with A (= 0.5), t in (= 100), and T w (= 10). Figures 3(a) and 3(b) 
show the time courses of p(t) and j(t,t) (= p(t,t)), respectively, with a set of parameters 
of a — 1, w — 0.5, P = 0.001 and N — 1, an input spike of I^(t) being shown at the 
bottom of Fig. 3(a). When an input spike is applied at t — 100, state variables of Xi(t) 
are randomized because independent noises have been added since t — 0. Solid and dashed 
curves in Fig. 3(a), which denote results of AMM and DSs, respectively, are practically 
identical. The dotted curve expressing p(t) of SDA is in fairly good agreement with that of 
DSs for r = 1, but the former completely disagrees with the latter for r > 2. We should note 
in Eqs. (20)-(23) that p(t) is decoupled from 7 (£, t) and p(t, t), then 7 (£, t) is independent of 
an external input I^ e \t). Figure 3(b) shows that time courses of 7 (t, t) of AMM and DS are 
almost identical for r = 0. For r = 1, result of AMM is underestimated compared to that 
of DS as discussed before. However, an agreement of the result of AMM with that of DS 
becomes better for r > 4. Results of 7 (£, t) of DS at large t (> 100) are in good agreement 
with the exact stationary solution of 7* shown in Fig. 2(c). 

N > 1 case 

We will discuss dynamics, in particular, the synchronization, of ensembles for N > 1. In 
order to monitor the synchronization, we consider the quantity given by [19] 



When all neurons are in the completely synchronous states, we get Xi(t) = X(t) for all i and 
then R(t) = 0. In contrast, in the asynchronous states, we get R(t) = 2(1 — 1/N)^{t,t) = 
Ro(t) because p(t,t) = j(t,t)/N [Eq. (16)]. The synchronization ratio S(t) is defined by 



which becomes 1 (0) for completely synchronous (asynchronous) states. The synchrony a s 
of the ensemble is defined by 



where the overline denotes the temporal average between times t± (=2000) and ti (=3000). 

Figures 4(a) and 4(b) show the time course of pit) and S(t), respectively, for various w 
with r = 10, P — 0.001 and iV = 10, solid and dashed curves denoting results of AMM and 
DS, respectively. For w — 0, pit) behaves as a simple relaxation process with the relaxation 
time of T r — 1/a — 1, while S(t) is vanishing. When a small, positive coupling of w = 0.4 is 
introduced, pit) shows the stair-like structure because of the positive delayed feedback. The 
synchronization ratio S(t) for w = 0.4 shows a gradual development as increasing t, but the 
magnitude of its averaged value of a s is very small (~ 0.015). When w is more increased 
to w — 0.8, the effective relaxation time for p(t) to return to the initial zero value becomes 
larger and a s becomes also larger. For w = 1.0, the effective relaxation time becomes infinity 



R (t) = ^2 E < [^(t) - x^t)} 2 >= 2[ 7 (M) -P(t,t)]. 



(26) 



[19] 




(27) 




(28) 
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and p(t) remains at the finite value of p(t) = 0.455. For w > 1, the divergence in p(t) is 
triggered by an input spike and S(t) tends to a fully synchronized state of a s = 1. On 
the contrary, for a small negative coupling of w = —0.4, p(t) shows an ostensibly quasi- 
oscillating state because of negative delayed feedback. With increasing the magnitude of 
negative w, the term showing this quasi-oscillation becomes longer. For w < —1.2, p(t) 
shows a divergent oscillation and S(t) tends to saturate at unity for t > 250. 

The w dependence of a s , which is the temporal average of S(t), is depicted in Fig. 5(a), 
where solid and dashed curves show results of AMM and DS with 1000 trials, respectively: 
error bars showing the root-mean-square (RMS) value of DS are within the radius of circles. 
Although a s is very small for | w |< 0.9, it is suddenly increased as | w | approaches the 
unity, where the divergence of the autonomous oscillation is induced as shown in Fig. 4(a). 

The delay time r plays an important role, as discussed before in the case of N = 1 [Fig. 
2(b)]. Figure 5(b) shows the r dependence of a s calculated by AMM (the solid curve) and 
DS with 1000 trials (the dashed curve) for w = 0.5, (3 = 0.001 and N = 10. We note that 
a s = 0.091 for r = is rapidly decreased with increasing r from zero, while it is almost 
constant at r > 4. This r dependence of a s resembles that of 7* for N = 1 shown in Fig. 
2(a). 

We have so far fixed the size of N, which is now changed. Figure 5(c) shows the N 
dependence of a s calculated in AMM (the solid curve) and DS (the dashed curve) for r = 10, 
w = 1.0 and (3 = 0.001. For N = 2, the synchronization of a s ~ 0.963 is nearly complete. 
With increasing N, however, a s is gradually decreased: a s — 0.824 and 0.340 for iV = 10 and 
100, respectively. 

Model calculations have shown that in linear Langevin ensembles with appropriate model 
parameters, an applied spike induces oscillations with divergent amplitudes. This is contrast 
with the nonlinear Langevin ensembles where stable oscillations with finite amplitudes are 
possible, as will be shown in the following Sec. IIIB. 



B. Nonlinear model 

Next we consider the nonlinear (NL) model in which F(x) and H(x) in Eq. (1) are given 

by 

F(x) = -ax, (29) 
H(x) = x- bx 3 . (a > 0, b > 0) (30) 

The NL model given by Eqs. (1), (29) and (30) with a = 0, b = 1, I^{t) = and iV = 1 
has been discussed in Ref. [10]. With the use of Eqs. (29) and (30), Eqs. (7)-(10) become 

M> = _ a ^ t ) + WUo{t _ T ) + /(e) (f) > (3i) 
at 

= -2a~f(t, t) + 2wu x (t - r)p(t, t - r) + (3 2 , (32) 



dt 

= _ 2ap (t, f) + 2w Ul (t - r)p(t, t-r) + ^, (33) 

dp(t, t — mr) 



dt 



= —2ap(t, t — mr) + wu\(t — (m + l)r)p(t, t — (m + l)r) 
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2 

+ wui(t - r)p(t -r,t- mr) + ^A(mr), (for m > 1) (34) 

with 

u (t) = p(t) - bp(tf - 3bp(t)~f(t, t), (35) 
Ul (t) = 1 - 3bp(t) 2 - 3&7(t, *). (36) 

For (3 = and I^ e \t) = 0, Eq. (31) has the stationary solution given by 

// = 0, for iy < a (37) 

for w > a (38) 

Linearizing Eq. (31) around p*, we get the condition for the stable stationary solution given 
by 

cos~ l (a/w) „ 
r < t c1 = = , for w < a (39) 

J (w 2 — a 2 ) 

cos~ 1 [a/(3a-2w)} . . 

< r c2 = , for w > 2a (40) 

^/[(3a - 2u>) 2 - a 2 ] 

Figure 6 shows the calculated w-t phase diagram of the NL model, showing the non- 
oscillating (NOSC) and oscillating states (OSC) with a = 1 and b = 1/6, which are adopted 
for a later comparison with the nonlinear model given by Eqs. (43) and (44) in Sec. IV. 
When an external spike given by Eq. (2) is applied to the NOSC state, the state is once 
excited and returns to the stationary state after the transient period, as will be discussed 
shortly [Figs. 8(a)-8(d)]. On the contrary, when a spike is applied to the OSC state, it 
induces the autonomous oscillation. 

Adopting parameters of w and r values shown by circles in Fig. 6, we have made 
calculations for the NL model by AMM and DS, whose results of /i(t) and S(t) are depicted 
in Figs. 7(a) and 7(b), respectively, with [3 = 0.001 and N = 10, solid and dashed curves 
expressing results of AMM and DS, respectively. Hereafter we show results of AMM6 whose 
validity for the NL model will be confirmed later [Fig. 9(b)]. Figure 7(a) shows that with 
increasing r, fi(t) shows the complicated time dependence due to delayed feedbacks. The 
time course of n(t) for iV = 10 is the same as that for N = 1 (results not shown). As was 
discussed in Sec IIIA, 7(t,t) and p(t,t) in the L model are independent of an input signal 
I^(t) because they are decoupled from pit) in Eq. (20)- (23). It is not the case in the NL 
model, where pit), j(t,t) and p(t,t) are coupled each other in Eqs. (31)-(34), and S(t) 
depends on /^(t). Figure 7(b) shows that, for example, in the case of r = 0, S(t) ~ 0.95 
at t ~ 100 is suddenly decreased to S(t) ~ at t = 100 by an applied spike, and then it is 
gradually increased to the stationary value of S* ~ 1.0 at t > 1000. This trend is realized 
in all the cases shown in Fig. 7(b). We note that an agreement of S(t) between AMM and 
DS is good for r = 0, 5 and 10, but not good for r = 1 and 2, just as in the case of the L 
model [Fig. 2(a)]. 

Figures 8(a)-8(d) show the time courses of p(t) and S(t) when the w value is changed 
along the horizontal dotted line in Fig. 6, solid and dashed curves denoting results of AMM 
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and DS, respectively. From comparisons among Figs. 4(a), 4(b), and 8(a)-8(d), we note that 
for | w | < 0.8, time courses of /j,(t) and S(t) of the NL model are similar to those of the L 
model. The difference between the L and NL models is, however, clearly realized in cases of 
| w |> 1.2. For w = —1.2, fi(t) in the NL model oscillates with the bounded magnitude [Fig. 
8(c)] while fj,(t) in the L model oscillates with divergent magnitude [Fig. 4(a)] although the 
oscillating period is the same (T = 22) for L and NL models. In contrast, S(t) in the NL 
model oscillates [Fig. 8(d)] while S(t) in the L model saturates at the unity [Fig. 4(b)]. For 
w = 1.2, fi(t) in the NL model starting from the stationary state with /z* = 1.0, is slightly 
modified by an input spike applied at t = 100 with a small magnitude of S* = 0.024 for S(t), 
whereas fj,(t) in the L model shows an unbounded oscillation and S(t) saturates at S* = 1. 
Figure 8(a) shows that as increasing w above 1.2, /z(t) shows a quasi-oscillation triggered by 
inputs, by which S(t) is increased at 110 ~ t ~ 130. For w > 2.0, the autonomous oscillation 
with a period of T = 22 is induced and S(t) is also oscillating with a period of T = 11. 

As discussed above, the oscillation is triggered by an input spike when parameters are 
appropriate. In order to study the transition between the NOSC and OSC states in more 
detail, we have calculated the quantity a defined by [31] 

^ = ^£[<^) 2 >-<^)> 2 ]> (41) 

i 

= - W) 2 + 7(M), (42) 

which becomes finite in the OSC state but vanishes (or is small) in the NOSC state, the 
overline denoting the temporal average [Eq. (28)]. 

Figure 9(a) shows a Q and a s calculated with changing w from 1.6 to 2.4 along the hor- 
izontal dotted line in Fig. 6 with r = 10, /3 — 0.001 and N = 10, solid and dashed curves 
denoting results of AMM6 and DS, respectively. The oscillation is triggered by an input 
spike when w exceeds the critical value of w c (~ 2.02). The transition is of the second 
order since a Q is continuously increased as (w — w c ) is increased. We should note that the 
peak in a s shows the fluctuation-induced enhancement at w ~ w c , which arises from an 
increase in the ratio of p(t, t)/~f(t, t) although both j(t, t) and p(t, t) are increased. When w 
exceeds about 2.3, the oscillation becomes irregular, which is expected to be a precursor of 
the chaotic state. 

For calculations of the NL model given by Eqs. (29) and (30), we have adopted AMM6, 
whose validity is examined in Fig. 9(b) showing a s for 1.8 < w < 2.6 when the level m is 
changed in AMMm: note that the result of m = 6 in Fig. 9(b) is nothing but the AMM 
result of a s in Fig. 9(a). The critical coupling for the NOSC-OSC transition calculated 
by AMMm for m = 1 — 3 is too large compared to that by DS (w c = 2.02) shown in Fig. 
9(a). For m = 4, we get a reasonable value of w c , but the peak value of a s at w = w c 
is too small. With furthermore increasing the m value, the w-dependence of a s becomes 
in better agreement with that of DS. The optimum value of m is expected to depend on 
the model parameters, the required accuracy and the ability of computer facility. Making a 
compromise among these factors, we have decided to adopt AMM6 in all our calculations. 
This choice of m = 6 has been confirmed by results of AMM which are in good agreement 
with those of DS [Fig 9(a)]. 

Figure 9(c) expresses the w dependence of a a and a s near the NOSC-OSC transition for 
a negative w when the w value is changed from -1.4 to -0.6 along the horizontal dotted line 
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in Fig. 6. The oscillation is triggered by an input spike when w is below the critical value 
of w c (~ -1.04). The fluctuation-induced enhancement is again realized in a s calculated by 
AMM and DS. 

The iV-dependence of a s has been studied for various w (~ 2) with r = 10 and (3 = 0.001, 
whose results are depicted in Fig. 10(a). It is noted that a s with w = 2.02 is larger than 
that with w = 2.04 for all N values, just as shown in Fig. 9(a). With increasing N, a s 
is gradually decreased for all w values. Although RMS values of a s in the NOSC state are 
small, they become considerable in the OSC states, which is due to oscillations in S(t) but 
not due to noises. In contrast, the iV-dependence of a is very small for the parameters 
investigated (results not shown). 

Figure 10(b) shows the (3 dependence of a Q and a s for w =2.0, 2.02 and 2.04 with r = 10 
and N = 10. For w = 2.02, which is just the critical coupling for (3 = 0.001 [Fig. 9(a)], a s 
is gradually decreased with increasing (3. For w = 2.00 (~ w c ), a s is little decreased with 
an increase in (3. In contrast, for w = 2.04, a s is first increased with increasing w and has a 
broad peak at (3 ~ f3 c where (3 C = 0.04 in AMM and 0.06 in DS. This broad peak in a s may 
suggest that the OSC state is suppressed by noises although a Q calculated in DS remains 
finite at (3 > (3 C , showing no signs for the NOSC-OSC transition. It may be possible that the 
emergence of the oscillation is not well represented by a Q defined by Eqs. (41) and (42) in 
the case of a large f3. The discrepancy between results of AMM and DS becomes significant 
with increasing (3, which is due to a limitation of AMM based on the weak-noise assumption. 

It has been shown that when model parameters of w, r, (3 and N are appropriate, 
the stable oscillations with finte magnitudes are induced in NL Langevin ensembles. The 
fluctuation-induced enhancement is realized in the synchrony near the second-order transi- 
tion between the NOSC and OSC states. 



which is referred to as the NL' model. Equations (43) and (44) were previously employed 
by Ikeda and Matsumoto [1] for a study on chaos in time-delayed systems. By using Eqs. 
(43) and (44), and the relation given by 



IV. CONCLUSIONS AND DISCUSSIONS 



We may adopt a nonlinear Langevin model given by Eq. (1) with 



F{x) = —ax, 
H(x) = sin(x), 



(a > 0) 



(43) 
(44) 



H {2n+1 \t) = (-l) n cos{x), 



(45) 
(46) 



we get DEs given by Eqs. (31)- (34) but with 




(48) 



(47) 
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where all contributions from n = to oo in Eqs. (11) and (14) are included. It is noted that 
Eq. (30) with b = 1/6 is an approximate form of Eq. (44) for a small x. Correspondingly, 
Uo(t) and U\{t) given by Eqs. (35) and (36) are approximate expressions of those given by 
Eqs. (47) and (48), respectively, for small p(t) and 7(£, £). 

Figure 11 shows the w dependence of a Q and a s of the NL' model. The NOSC-OSC 
transition occurs at w c = 2.29 above which a is continuously increased and where a s has a 
peak. Although the induced oscillation is regular for 2.29 < w ~ 2.64, it becomes irregular 
for w ~ 2.64, which may lead to the bifurcation and chaos [1]. We note from Fig. 9(a) 
and 11 that the w dependence of the NL model given by Eqs. (29) and (30) is similar 
to that of the NL' model given by Eqs. (43) and (44) although the critical coupling for 
the NOSC-OSC is different between the two models. The w — r phase diagram for NL' 
Langevin model is almost the same as that for NL Langevin ensembles shown in Fig. 6. 
Recently the phase diagram has been experimentally obtained for a coupled pair of the 
Plasmodium of the slime mold, physarum polycephalum, where the coupling strength and 
delay time are systematically controlled [32]. The observed phase diagram is not dissimilar 
to our w — r phase diagrams for NL and NL' Langevin models as far as unentrained and 
in-phase oscillating states are concerned. 

Quite recently, Huber and Tsimring (HT) [17] have discussed an alternative nonlinear 
Langevin ensembles given by Eq. (1) with 

F(x) = x-x 3 , (49) 
H(x) = x, (50) 

for l( e \t) = 0, which expresses interconnected bistable systems with delays [14]. By using 
DSs and analytical methods based on the Gaussian and dichotomous approximations, HT 
have discussed the coherence resonance and multistability of the system. When we apply 
our approximation to this nonlinear model, Eqs. (7)-(10) become 

3/i(t) 7 (t, t) + wp(t - t) + I^it), (51) 

37(t, i)] 7 (i, t) + 2wp(t, t-r) + P\ (52) 

3 7 (t, t)]p(t, t) + 2wp(t, t-r) + ^, (53) 

mr)]p(t, t — mr) + wp(t, t — (m+ l)r) 

+ p(t-r,t- mr) + ^A(mr), (for m > 1) (54) 

where gi(t) = l—3p(t) 2 — 3 , ~f(t,t). In their Gaussian approximation, HT have employed Eqs. 
(51) and (52) with p(t,t — r) = 0, discarding Eqs. (53) and (54). It has been claimed that 
the Gaussian approximation is not adequate near the transition point between the ordered 
and disordered states, although a dichotomous theory yields a fairly good description. This 
might be due to a neglect of the higher-order contributions in Eqs. (52)- (54), which are 
expected to play important roles, in particular, near the transition point, as our calculations 
have shown [Fig. 9(b)]. 



dp(t) 

dt 
drf(t,t) 

dt 
dp(t,t) 

dt 

dp(t, t — mr) 
Jt 



p(t)-p(tf 
2[1 - 3/i(t) 2 
2[1 -3/i(t) 2 

[gi(t)+gi(t 
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In summary, we have proposed a semi-analytical approach for a study of dynamics of 
stochastic ensembles described by linear and nonliner Langevin models with delays. Advan- 
tages of our method are (a) the synchronization in ensembles may be discussed by taking 
into account correlations of local and global variables, (b) the recursive DEs terminated at 
finite m (~ 6) yield fairly good results compared to those of DSs, and (c) our method is free 
from the magnitude of time delays though the noise intensity is assumed to be weak, which 
is complementary to SDA [9]. The proposed method is expected to be useful not only to 
Langevin ensembles but also to more general stochastic ensembles with delays. Although 
our method is applicable to the system with an arbitrary size of N, it is better applied to 
larger system because of its mean-field nature. It should be noted that the number of DEs 
to be solved for A-unit stochastic Langevin model is NN r in DS with iV r trials, while it is 
(m + 3) in AMMm. The ratio between the two numbers becomes NN r / (m + 3) ~ 1000, for 
example, for N = N r = 100 and m — 6. Actually this reflects on the ratio of the speed for 
numerical computations by using the two methods. Taking these advantages of our method, 
we have applied it to ensembles described by FN neuron model with delayed couplings to 
study their dynamics and synchronization, which are reported in a following paper [26]. 
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APPENDIX A: DERIVATION OF EQS. (7)-(10) 

Assuming that the noise intensity (5 is small, we express Eq. (1) in a Taylor expansion 
of Sxi as 

= t + ^1 - r Y + m + /'•'(«). (Al) 

at e=Q i. iv . e=0 i. 

where F®(t) = F®(n(t)) and H^(t) = H^(fi(t)). Equations (3), (4) and (Al) yield DE 
for means of //(£) as 

+ ^EEi; £ ^%^<^(t-r)'>+/MW. (A2) 

When we adopt the Gaussian decoupling approximation, averages higher than the second- 
order moments in Eq. (A2) may be expressed in terms of the second-order moments given 
by 

< 5x±, ...,5xg > = n fem < 5xkSx m >, for even £, 

all parings 

= 0, for odd £, (A3) 
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where the summation is performed for all (£ — 1)(£ — 3). ...3 ■ 1 combinations. With the use 
of the Gaussian decoupling approximation given by Eq. (A3), Eq. (A2) becomes 

d»(t)_ 1 V f^ (2n) W 5 <Sx(t?> n 

+ £ E E E g( ^n)7 r) g2w < ^ ~ r)2 >n +/(e)(t) ' (A4) 

where -B 2n = (2n — l)(2n — 3) • -3 • 1. Adopting the mean-field approximation given by 

< 5x, t (t) 2 > n ~ 7 (t, t)™" 1 < 5^(t) 2 >, (A5) 

we get 

dfi(t) ~ F( 2 ")(t) fi(t,t)\ n 



dt 



t ^ (2M)" + . | (1(^1))" + „ (t) , (A6) 



which yields Eq. (7), (11) and (13). 

From Eqs. (Al) and (A4), we gt DEs for dSxi(t)/dt as 

dSxi(t) dxiit) dji{t) 
dt dt dt 

~to (t) (2n+l)! + „V W V (2n)! ~" 2™ „! , 

w - r) 2 " +1 



(A7) 



By using Eqs. (5), (A3) and (A8), we get DEs for dj(t,t)/dt as 

^4 ?<Mi) ^> (A9) 

+ ^ZSE ^+'iV' * Mi)M ' " T)2 " +1 > * fa,(%(t) " 



00 /a - ^ Y+ _ ^2ra 



^(t-r) 2n 7(t-r,t-r) , 



+ |S E E E H ( ^^J ) B 2n+2 < 5 Xl {t)5x 3 {t - r) >< 6 Xj (t - rf > 



+ /3 2 . (A10) 
With the use of the mean-field approximation given by Eq. (A5), Eq. (A10) reduces to 
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*y(M) -2 7 (t,of; i?<an+1)( * ) ^ {t,t)Y 



dt ' ^ n\ 



+ 2mp(M _ T) f^^)(2(^l))" + ^ (A11 ) 

leading to Eq. (8), (12) and (14). Calculations of dp(t, t)/dt and dp(t, t—mr)jdt are similarly 
performed by using the relation: 

dp(t,t — mr) 1 v^v-^ . ..dxAt — mr) dxAt) , . /»,^s 

* 3 

In the process of calculating dp(t,t — mr)/dt, we get new correlation functions given by 

S(t, t - mr) = — ^ < fe,(t)&(t - mr) >, (A13) 

i 

S(t -mT,t) = —^2 < Sxi(t - mT)ti(t) > . (A14) 



By using the method of steps in Ref. [13], we get 



S(t, t-mr) = S(t - mr, t) = (^j A(mr), (A15) 
which leads to Eq. (10). 

APPENDIX B: THE SMALL-DELAY APPROXIMATION 

We apply the small-delay approximation (SDA) first proposed in Ref. [9] to our model 
given by Eqs. (1) and (2). When r is small, we may expand Eq.(l) for iV = 1 as x(t — r) ~ 
x(t) — r dx(t)/dt to get 

^ ~ F(x(t)) + w (H(x(t)) - r H '(x(t))^j + ^{t) + &\t). (Bl) 

Using Eq. (Bl), we get DEs for pit) and j(t,t) given by 

^ = [1 - wrh^lgoit) + wu (t) + I {e \t)}, (B2) 

" ; " ' ' = 2[1 - wrh^Wg^t) + wu 1 {t)} 1 {t,t) + [1 - wrh^t)} 2 ^, (B3) 



dt 

where h^t) = H'{p{t)). For the L model given by Eqs. (1), (2), (18) and (19), Eqs. (B2) 
and (B3) become 

- (l-«;r)[(-o + u;Ht) + / (e) (i)] 1 (B4) 



dt 

dl(t, t) _ _ wr ^_ a + w ^ t ) + (1 _ W1 f(p. (B5) 



dt 

The r dependence of the stationary solution of 7* is shown by the dotted curve in Fig. 2(b) 
The time course of pit) is plotted by doted curves in Fig. 3(a). 
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FIGURES 



FIG. 1. The recursive structure of equations of motions p(t,t) in AMM, arrows denoting the 
mutual dependence (see text). 

FIG. 2. (a) The stationary solution of "f(t, t), 7*, of the linear (L) model given by Eqs. (18) and 
(19) calculated in AMMm with changing the level m (solid curves) and in the exact calculation 
(dashed lines) for various r with a = 1, w = 0.5, (3 = 0.001 and N = 1. (b) The r dependence of 7* 
in AMM6 (the solid curve), SDA (the dotted curve) and exact calculations (the dashed curve) for 
a = 1, w = 0.5, (3 = 0.001 and N = 1. (c) The w dependence of 7* for various r in AMM6 (solid 
curves) and exact calculations (dashed curves) for [3 = 0.001 and N = 1. Results are multiplied by 
a factor of 10 6 , and those of r = 5, 2, 1, and in (a) and (c) are successively shifted upwards by 1. 

FIG. 3. (color online). The time course of (a) fi(t) and (b) "y(t,t) of the L model for an 
applied single spike shown at the bottom frame in (a), calculated in AMM (solid curves), the 
small-delay approximation (SDA; dotted curves) and direct simulations (DS; dashed curves) with 
a = 1, w = 0.5, (3 = 0.001, and N = 1: results of fi(t) for AMM and DS are indistinguishable, and 
results of SDA are shown only for r < 2 (see the text). 

FIG. 4. (color online). Time courses of (a) /i(t) and (b) S(t) of the L model calculated by AMM 
(solid curves) and DS (dashed curves) for various w with r = 10, (3 = 0.001 and N = 10: results 
of (i(t) for AMM and DS are indistinguishable. The chain curve at the bottom of (a) expresses an 
applied input spike. 

FIG. 5. (a) The w dependence of a s , the temporal average of S(t), of the L model calculated 
in AMM (the solid curve) and DS with N r = 1000 (the dashed curve) for r = 10, (3 = 0.001 and 
N = 10. (b) The r dependence of a s for w = 0.5, (3 = 0.001 and N = 10 calculated by AMM (the 
solid curve) and DS with iV r = 1000 (the dashed curve), (c) The dependence of a s for w = 0.5, 
r = 10 and (3 = 0.001 calculated by AMM (the solid curve) and DS (the dashed curve); N r = 100 
for N = 50 and 100, and N r = 1000 otherwise. Error bars denote RMS values of DS. 

FIG. 6. The w-t phase diagram of the nonlinear (NL) model given by Eqs. (29) and (30) with 
a = 1, g = 1/6 and (3 = 0, showing the non-oscillating (NOSC) and oscillating (OSC) states. 
Calculations whose results are depicted in Fig. 7, are performed for sets of parameters shown by 
circles. Along the horizontal dotted lines, the w value is changed for calculations shown in Figs. 8 
and 9(a)-9(c) (see text). 

FIG. 7. (color online). The time course of (a) /i(i) and (b) S{t) of the NL model for an 
applied single spike shown at the bottom frame in (a), calculated in AMM (solid curves) and DS 
(dashed curves) for a = 1, w = 1.0, (3 = 0.001, and N = 10: results of fi(t) for AMM and DS are 
indistinguishable. 
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FIG. 8. (color online). Time courses of (a) fj,(t) and (b) S(t) for 1.2 < w < 2.3, and (c) ji{t) and 
(d) S(t) for -1.2 < w < 1.0, of the NL model calculated by AMM (solid curves) and DS (dashed 
curves) with r = 10, f3 = 0.001 and N = 10: results of fj,(t) for AMM and DS are indistinguishable. 
Chain curves at bottoms of (a) and (c) express applied input spikes. 

FIG. 9. (a) The w dependence of a Q and a s for 1.6 < w < 2.4 of the NL model, (b) The 
w dependence of a s for 1.8 < w < 2.6 with different level m in AMMm (see text), (c) The w 
dependence of a and a s for —1.4 < w < —0.6. Solid and dashed curves in (a) and (c) express 
results of AMM6 and DS calculated with r = 10, f3 = 0.001 and N = 10. Note that the result with 
m = 6 in (b) corresponds to the AMM result of a s in (a). Errors bars expressing RMS values are 
not shown for a clarity of figures (see Fig. 10) 

FIG. 10. (a) The N dependence of a s for r = 10 and (3 = 0.001, and (b) the (3 dependence of 
a s for r = 10 and N = 10 of the NL model, calculated by AMM (solid curves) and DS (dashed 
curves) with w = 2.04 (triangles), 2.02 (circles), 2.0 (circles), 1.95 (diamonds) and 1.90 (inverted 
triangles). Errors bars expressing RMS values of a s of DS, are significant in the OSC state because 
of the oscillation of S(t). 

FIG. 11. The w dependence of a and a s of the NL' model given by Eqs. (43) and (44), 
calculated by AMM (solid curves) and DS (dashed curves) with r = 10, (5 = 0.001 and TV = 10 
(see text). Errors bars expressing RMS values are not shown for a clarity of the figure. 
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